Enteropathogen seroepidemiology among children in low-resource settings

Asembo, Kenya antibody profiles, Figure 5

Notebook Summary

This script plots individual antibody trajectories in the Asembo, Kenya cohort from enrollment to follow-up. 199 children were measured at both timepoints.

Script preamble

#-----------------------------
# preamble
#-----------------------------
library(here)

## here() starts at /Users/benarnold/enterics-seroepi

here()

## [1] "/Users/benarnold/enterics-seroepi"

# load packages
library(tidyverse)

## ── Attaching packages ─────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 3.0.0     ✔ purrr   0.2.4
## ✔ tibble  1.4.2     ✔ dplyr   0.7.4
## ✔ tidyr   0.8.0     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0

## ── Conflicts ────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

# bright color blind palette:  https://personal.sron.nl/~pault/ 
cblack <- "#000004FF"
cblue <- "#3366AA"
cteal <- "#11AA99"
cgreen <- "#66AA55"
cchartr <- "#CCCC55"
cmagent <- "#992288"
cred <- "#EE3333"
corange <- "#EEA722"
cyellow <- "#FFEE33"
cgrey <- "#777777"

#-----------------------------
# load the formatted data
# created with 
# asembo-enteric-ab-data-format.Rmd -->
# asembo-enteric-ab-distributions.Rmd
#-----------------------------
dl <- readRDS(here("data","asembo_analysis2.rds"))

Identify incident changes

Identify children whose status changed between enrollment and follow-up. Those who changed from negative to positive are seroconverters (seroi below), and those who changed from positive to negative are seroreverters (seror below).

#-----------------------------
# identify incident 
# seroconversions and reversions
# based on crossing the
# seropositivity cutoff
# and
# based on a 4-fold change in MFI
# to above the cutoff (seroconversion)
# or starting above cutoff (seroreversion)
#-----------------------------
dlp <- dl %>%
  group_by(antigen,antigenf,childid) %>%
  mutate(seroposA=ifelse(time=="A",seropos,NA),
         seroposA=max(seroposA,na.rm=TRUE),
         seroposB=ifelse(time=="B",seropos,NA),
         seroposB=max(seroposB,na.rm=TRUE),
         seroi=ifelse(seroposB==1 & seroposA==0,1,0),
         seror=ifelse(seroposB==0 & seroposA==1,1,0),
         
         mfiA=ifelse(time=="A",logmfi,NA),
         mfiA=max(mfiA,na.rm=TRUE),
         mfiB=ifelse(time=="B",logmfi,NA),
         mfiB=max(mfiB,na.rm=TRUE),
         seroi4fold=ifelse(logmfidiff>log10(4) & mfiB>serocut,1,0),
         seror4fold=ifelse(logmfidiff< - log10(4) & mfiA>serocut,1,0)
  ) 

Antibody MFI trajectories

Longitudinal antibody trajectories of individual children for all antigens measured, colored by the change in antibody levels (increases and decreases)

#-----------------------------
# Plot individual antibody
# trajectories, colored
# by increases and decreases
#-----------------------------

# create an change category factor for plot aesthetics
dlp <- dlp %>%
  mutate(mficat=factor(ifelse(logmfidiff>0,"Increase","Decrease"),levels=c("Increase","Decrease")))

# convert time to numeric (for plotting)
dlp <- dlp %>%
  mutate(svy=ifelse(time=="A",0,1))

log10labs <- c( 
  expression(10^0),
  expression(10^1),
  expression(10^2),
  expression(10^3),
  expression(10^4)
)

ppq <- ggplot(data=filter(dlp,!is.na(mficat)),aes(x=svy,y=logmfi,group=factor(childid),color=mficat)) +
  # geom_point(alpha=0.2) +
  geom_line(alpha=0.2) +
  # geom_hline(aes(yintercept=mixcut),linetype="dashed") +
  
  scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
  scale_y_continuous(breaks=0:4,labels=log10labs) +
  scale_color_manual(values=c(corange,cteal,"gray80"),
                     guide=guide_legend(title="MFI change:",override.aes = list(alpha=1))
                     ) +
  facet_wrap(~antigenf,nrow=6,ncol=2) +
  labs(x="Measurement",y="Luminex response (MFI-bg)") +
  theme_minimal() +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.ticks.y=element_line(color="gray40"),
    legend.position = "top"
    )

ppq

Create the same figure of individual antibody trajectories between measurements, but limit it to antigens with seroincidence measures, and color the trajectories by seroconversion/reversion status.

#-----------------------------
# pull the incidence info
# and just merge it back
# to the full longidudinal data
#-----------------------------
dl2 <- dlp

# create an incidence category factor for plot aesthetics
# for seroconversion based on seropositivity only
dl2$icat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$icat[dl2$seroi==1] <- "Seroconversion"
dl2$icat[dl2$seror==1] <- "Seroreversion"
dl2$icat[dl2$seroi==0 & dl2$seror==0] <- "No change"


# create an incidence category factor for plot aesthetics
# for seroconversion based on 4-fold changes in MFI
dl2$i4cat <- factor(NA,levels=c("Seroconversion","Seroreversion","No change"))
dl2$i4cat[dl2$seroi4fold==1] <- "Seroconversion"
dl2$i4cat[dl2$seror4fold==1] <- "Seroreversion"
dl2$i4cat[dl2$seroi4fold==0 & dl2$seror==0] <- "No change"

# create an incidence category factor for plot aesthetics
# for seroconversion comparing both methods
dl2$icat_comp <- factor(NA,levels=c(">4-fold increase, across cutoff",">4-fold decrease",">4-fold increase, above cutoff","<4-fold change"))
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==1] <- ">4-fold increase, across cutoff"
dl2$icat_comp[dl2$seror4fold==1] <- ">4-fold decrease"
dl2$icat_comp[dl2$seroi4fold==1 & dl2$seroi==0] <- ">4-fold increase, above cutoff"
dl2$icat_comp[is.na(dl2$icat_comp)] <- "<4-fold change"


pp <- ggplot(data=filter(dl2,!is.na(icat)),aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
  # geom_point(alpha=0.2) +
  geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
  geom_line(alpha=0.2) +
  
  scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
  scale_y_continuous(breaks=0:4,labels=log10labs) +
  scale_color_manual(values=c(corange,cteal,cmagent,"gray80"),
                     guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
                     ) +
  facet_wrap(~antigenf,nrow=6,ncol=2) +
  labs(x="Measurement",y="Luminex response (MFI-bg)") +
  theme_minimal() +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.ticks.y=element_line(color="gray40"),
    legend.position = "top" 
    )

pp

Figure 5, antibody trajectories

#-----------------------------
# final figure, excluding cholera
# due to cross-reactivity with ETEC
#-----------------------------
relev <- levels(dl2$antigenf)[c(1:6,11,7,9,10)]
dl3 <- dl2 %>%
  filter(!is.na(icat) & antigen!="cholera") %>%
  ungroup() %>%
  mutate(antigenf=factor(antigenf,levels=relev))


pp2 <- ggplot(data=dl3,aes(x=svy,y=logmfi,group=factor(childid),color=icat_comp)) +
  # geom_point(alpha=0.2) +
  geom_hline(aes(yintercept=serocut),linetype="dashed",size=0.3) +
  geom_line(alpha=0.2) +
  
  scale_x_continuous(limits=c(-0.25,1.25),breaks=c(0,1),labels=c("Enrollment","Follow-up"))+
  scale_y_continuous(breaks=0:4,labels=log10labs) +
  scale_color_manual(values=c(corange,cteal,cmagent,"gray80"),
                     guide=guide_legend(title="MFI change:", override.aes = list(alpha=1), ncol=2,nrow=2),
                     ) +
  facet_wrap(~antigenf,nrow=6,ncol=2) +
  labs(x="Measurement",y="Luminex response (MFI-bg)") +
  theme_minimal() +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.ticks.y=element_line(color="gray40"),
    legend.position = "top" 
    )

pp2

# save PDF and TIFF versions
ggsave(here("figs","Fig5-asembo-ab-trajectories.pdf"),plot=pp2,device=cairo_pdf,width=6,height=10)
ggsave(here("figs","Fig5-asembo-ab-trajectories.TIFF"),plot=pp2,device="tiff",width=6,height=10)

Session Info

sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2  forcats_0.3.0   stringr_1.3.1   dplyr_0.7.4    
##  [5] purrr_0.2.4     readr_1.1.1     tidyr_0.8.0     tibble_1.4.2   
##  [9] ggplot2_3.0.0   tidyverse_1.2.1 here_0.1       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.18     cellranger_1.1.0 pillar_1.2.2     compiler_3.5.0  
##  [5] plyr_1.8.4       bindr_0.1.1      tools_3.5.0      digest_0.6.16   
##  [9] lubridate_1.7.4  jsonlite_1.5     evaluate_0.10.1  nlme_3.1-137    
## [13] gtable_0.2.0     lattice_0.20-35  pkgconfig_2.0.2  rlang_0.2.2     
## [17] psych_1.8.4      cli_1.0.0        rstudioapi_0.7   yaml_2.1.19     
## [21] parallel_3.5.0   haven_1.1.1      withr_2.1.2      xml2_1.2.0      
## [25] httr_1.3.1       knitr_1.20       hms_0.4.2        rprojroot_1.3-2 
## [29] grid_3.5.0       glue_1.2.0       R6_2.2.2         readxl_1.1.0    
## [33] foreign_0.8-70   rmarkdown_1.9    modelr_0.1.2     reshape2_1.4.3  
## [37] magrittr_1.5     backports_1.1.2  scales_1.0.0     htmltools_0.3.6 
## [41] rvest_0.3.2      assertthat_0.2.0 mnormt_1.5-5     colorspace_1.3-2
## [45] stringi_1.2.2    lazyeval_0.2.1   munsell_0.5.0    broom_0.4.4     
## [49] crayon_1.3.4